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A FINITE ELEMENT APPROACH TO THE STRUCTURAL 
INSTABILITY OF BEAM COLUMNS, 

FRAMES, AND ARCHES 


SUMMARY 


From the principle of virtual displacements and the bifurcation theory of 
elastic stability a stiffness matrix is developed for a beam column element with 
shear, moment, and axial load applied to the ends (nodes) of the element and a 
uniformly distributed load applied along the span of the element. The stiffness 
matrix developed relates the forces applied at the nodes to the displacements of 
the nodes and is a function of the magnitude of the applied load. Three types of 
behavior of the applied loads are considered as separate cases. These are: the 
loads remain normal to the deformed element; the loads remain parallel to their 
original direction; and, the loads remain directed toward a fixed point. 

A method is given for using the element stiffness matrix to predict the 
buckling load for a structure which may be represented by beam column elements. 
As an example, the buckling load of an arch for each of the three load-behavior 
cases is calculated and compared to known solutions. 


INTRODUCTION 


In recent years interest in the solution of structural problems by the 
finite element method has greatly increased. This interest can be attributed to 
the fact that the digital computer has made possible the solution of a great num- 
ber of complex, yet practical, problems by this method. 

Those active in publishing results in the field include a group at the 
University of Washington and The Boeing Company [ 1—9] , J. H. Argyris at the 
University of London [10, 11] , and R. H. Gallagher and his associates at Bell 
Aerosystems Company [12-14]. 

Basically the method consists of representing an actual structure by an 
idealized structure made up of elements for which relationships between the 
forces applied to points (nodes) on the boundary of the element and the displace- 
ments of the nodes is known. This relationship is conveniently represented by 


an element stiffness matrix. If all of the element stiffness matrices are trans- 
formed to a common (global) coordinate system, the sum of forces applied to 
the elements meeting at each node is equated to the externally applied force at 
that node, and the displacements of elements which meet at each node are equated, 
a set of equations relating the externally applied forces to the nodal displacements 
results. When this set is expressed in matrix form, the matrix relating the ex- 
ternally applied forces to the nodal displacements is called the master stiffness 
matrix. Then for any given combination of known applied external forces and 
nodal displacements the above relationships may be algebraically manipulated to 
determine internal loads and displacements throughout the structure. 

Once a sufficient variety and quality of element stiffness matrices are 
available to idealize a particular structure, the remainder of the above procedure 
to solve a particular problem becomes tedious but quite mechanical. Hence the 
advantage of the computer in solving such problems is obvious. 

Although the use of such techniques in solving complex structural prob- 
lems is now an everyday occurrence, a great deal of developmental activity in 
the field still persists for several reasons. First, there is no such thing as a 
correct stiffness matrix for a particular element which excludes other stiffness 
matrices for that same element as incorrect. A number of correct stiffness 
matrices may be developed for the same element if different assumptions are 
made for the derivation. Each of these matrices may have its own advantages 
for some particular application and all of them may converge to identical answers 
if the structure to be analyzed is broken into small enough elements. It is the 
search for the most advantageous element for particular applications or for all- 
inclusive elements which gives rise to much of the current research. 

Second, many derivations have been based on geometrical considerations 
only and have sometimes led to incorrect stiffness matrices or matrices with 
poor convergence characteristics. The recent trend has been to base derivations 
on basic elasticity principles and to use more care in assuring that adjacent 
element deformations conform to each other not only at the nodes but along the 
complete boundary [5, 15-19]. 

Third, extension of the finite element techniques to include a larger class 
of problems has received a great deal of recent attention. In particular, the 
solution of nonlinear problems including large deformations and structural 
instability has been of interest [1-3, 5, 8, 9, 13-15, 17, 20]. 
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This report was motivated to some extent by all three of the above con- 
siderations. The objective was to develop a beam column element, beginning 
with fundamental principles, which would account for the applied load behavior 
and could be used for structural stability analysis. 

It is assumed that the reader is familiar with the conventional matrix and 
index notation used throughout the report. Since there is no stringent space 
limitation, the work is presented in an unusual amount of detail with the pre- 
sumption that the development could serve as a guide for the development of 
curved or three-dimensional elements. 


GENERAL THEORY 
Basic Assumptions and Limitations 


The usual Euler-Bernoulli beam assumptions are made for developing 
the beam column element stiffness matrix. These are basic engineering assump- 
tions and may be listed as follows: 

1. The material of the element is homogeneous and isotropic. 

2. Plane sections remain plane after bending. 

3. The stress-strain curve is identical in tension and compression. 

4. Hooke's law holds. 

5. The effect of transverse shear is negligible. 

6. The deflections are small compared to the cross sectional, dimensions. 

7. The loads act in a single plane passing through a principal axis of 
inertia of the cross section. 

8. No initial curvature of the element exists. 

9. No local type of instability will occur within the element. 

10. Loads are applied quasi-statically. 
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The above assumptions limit the application of results to structures that 
lie in a plane and are loaded in the same plane. Out-of-plane and torsional 
buckling are not treated. 


Method of Analysis 


The work done during a virtual displacement of the element is equated 
to zero, as presented by Hoff in Reference 21 and derived more rigorously by 
Langhaar in Reference 22, to obtain the equilibrium equations for the element. 
The bifurcation concept of elastic stability as presented by Novozhilov in Refer- 
ence 23, Chapter V, is used to postulate that two possible sets of displacements 
which satisfy the equilibrium equations may exist under the same magnitude of 
external load if the magnitude is such that a structure is unstable. Each of these 
sets of displacements is substituted in turn into the equilibrium equations. The 
resulting sets of equations are combined to obtain a relationship between the 
nodal forces and the nodal displacements during buckling. When placed in matrix 
form this relationship becomes 


{Q 1 } = [[K] + [G] + [L]] {q 1 } = {0} 


( 1 ) 


where Q 1 is the column matrix of nodal forces , q 1 is the column matrix of nodal 
displacements, K is the usual beam element stiffness matrix, G is the geometric 
stiffness matrix, and L is the load behavior stiffness matrix. G and L both 
contain the magnitude of external load as a factor. Thus an element stiffness 
matrix can be derived which is a function of the applied load including the applied 
load behavior. 

A number of such elements may be combined to represent a particular 
structure, so that equation ( 1) applies to the entire structure and K + G + L is 

n^J 

the master stiffness matrix. Boundary conditions may be applied to reduce the 
size of the master stiffness matrix. For instability to exist the determinant of 
the master stiffness matrix must vanish. Hence, an eigenvalue problem is 
formulated where the eigenvalues are the magnitudes of the applied load at which 
the structure is unstable. 
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Development of Element Stiffness Matrix 


Description of Element. Consider a beam column element as shown in 
Figure 1 which is subjected to nodal forces and moments and a uniformly dis- 
tributed load, p, applied along its length. It is desired to determine a stiffness 
matrix for the element which may be used to calculate the stability of structures 
made up of such elements. The stiffness matrix is to account for the fact that 
the components of the nodal forces or distributed load, or both, may be a function 
of the element displacements. Three cases are to be considered: the loads 
remain normal to the deformed element; the loads remain parallel to their 
original direction; and, the loads remain directed toward a fixed point. 



FIGURE 1. BEAM COLUMN ELEMENT 

Displacement Functions . If it is assumed that the lateral displacement 
of the beam element of Figure 1 may be represented by 

w = a j + a 2 x + a 3 x 2 + a 4 x 3 (2) 

and the longitudinal displacement is given by 

u = a 5 + Q! e x , ( 3) 
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then the six unknown constants, c^, may be determined in terms of the nodal 
displacements from the element boundary conditions , 


w 

= Wj 

at 

x = 0 

w ’x 

= -e t 

at 

x = 0 

w 

= w 2 

at 

x = L 

W, 

X 

= -H 

at 

x = L 

u 

= Ui 

at 

X 

II 

o 

u 

= u 2 

at 

x = L 


The resulting displacement functions are 



(4) 


= W I <Pi(x) + W 2 0 2 ( x ) + 03 (x) + $2 04 (x) 


= w. 0.(x) ( 5) 

1 1 

= Uj yj(x) + u 2 y 2 (x) 

* u. y^x) . (6) 
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These displacement functions have been used by others for the beam 
column element (see Reference 5, for example) . However, it should be men- 
tioned that an inconsistency arises when the cubic function is used for the lateral 
displacement of the present element which has a distributed load. This is illus- 
trated by taking the third derivative of w, which should be the equation for the 
shear load on the element, and observing that a constant results. But the ele- 
ment has a distributed load, and the shear should obviously be a linear function, 
not a constant. Since the change in shear over the length of an element becomes 
negligible in the limit as the element becomes smaller and smaller, this incon- 
sistency is not unacceptable and it will be seen that adequate results are obtained. 
The correct fourth order function could not be used for w because there is no 
oilier boundary condition available for evaluating another constant in equation ( 2) . 
This will be discussed further in the conclusions. 

Development of Equilibrium Equations . The element equilibrium equa- 
tions will now be developed from the principle of virtual displacements, 


SU - 6V = 0 (7) 

where SU is the change in strain energy and <5V is the work done by the external 
forces during a virtual displacement. 

For the uniaxial state of stress assumed to exist in a beam the change 
in strain energy is given by 


SU = f a 6e dv (8) 

J X X 
V 

Here, e is the strain of the beam at any point of the cross section and is given 
by x 


e = e + z k (9) 

x xx xx 


where e is the strain of the beam mid-surface, z is the distance of the point in 
xx 

the beam from the mid-surface, and 


k = w, 

XX XX 


( 10 ) 
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is the beam curvature. The nonlinear strain-displacement relationship 


e 

xx 





(ID 


is used for the mid- surface strain. The u, 2 term in this expression is not 

known to have been used in previous derivations of nonlinear beam element 
stiffness matrices. Substituting equations (9), (10), and (11) into equation 
( 8 ), 


<5U = f / Ee dxdA = f f E(e + ?k ) (6e + z&K )dxdA 

J A J L X X J ^ XX XX XX XX 


( 12 ) 


«u = / / 

A L 


Ee <5 e + Eze 5k + Ezk <5e + Ez 2 k 5k dxdA 

xx xx xx xx xx xx xx xx 


J zdA = 0 , J z 2 dA * I 
A A 


6U=/ 

L 


AEe <5e + El k 5 k 

xx xx xx xx 


dx 


C13) 


where 


Se = <5u, + 6 

xx x 


5 k = <5w, 

xx xx 


(l w 'x) + 6 (| u ’x) * 5u - 


+ w, <5w, + u, 6u, 

X X X X X 


5U “ L EA ( U ’ X + 2 W 'x + 2 U 'x) ( <5 % + w> x aw> x + U 'x Su 'x) +EIw 'xx 6w 'xx] dx 



dx 


<5U = 




+ EA 


( U -X w 'x) 


<5w, 


+ EIw , <5w , 

xx xx 


+ higher order terms. (14) 

The virtual work of the external forces acting on the element is given by: 


6V = F <5u. + F <5w. + M$0. + f p 6w + p Su dx (15) 

XllZllll£|_*Z X J 

where the force components depend upon the load behavior. Three load behavior 
cases will be considered. 


Case I. Loads Remain Normal to the Deformed Element ( Fig. 2) 



F . 

= F . 

+ F . 

e. 

XI 

XI 

Zl 


F . 

= F . 

- F . 

9. 

Zl 

Zl 

XI 

l 


M. = M. 

l l 


(16) 


Then 


<5V 


( F . + F .0. 

XI ZX 1 


6u. + 

l 



- f .e. 


XI 1 



+ M.60. 

l l 


+ 


/ pdw - pw, 

L L X 



dx 


(17) 
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FIGURE 3. LOADS REMAIN PARALLEL TO ORIGINAL DIRECTION 


6V = F . Su. + F . 6w. + M. <50. + / p<5wdx ( 19 ) 

xx x zx x x x 

Case III. Loads Remain Directed Toward a Fixed Point ( Fig. 4) 


F . 

R 

= F . + 

F . 

u i 

XI 

XI 

Zl 

R 

F . 

= F . - 

F . 

Wj 

Zl 

zl 

XI 

1 


M. = M. 
x x 


( 20 ) 
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Substituting equations ( 14) , (17) , ( 19) , and (21) into equation (7) , we 
obtain the following results: 

Case I. 

/ [ EA (% + h-l + 1 1 a -x) ,5u - x + EA ( u -x w ' X ) aw - x 

dx = [ F . + F . G. Uu; +(f . - F .6. W. 

4J y xi zi _i / i l zi xi jl J 1 


+ EIw, 6w, 

xx xx 


+ M.<50. + J p<5w - pw, x <5u, 
L 


dx 


( 22 ) 


Case II. 


/[ EA (% + I W 'x + f U 'x ) 6u 'x + EA (“'x W -x) 5w 'x 

+ EIw, 6w, dx = F . 6u. + F . <5w. + M.60.+ f pdwdx 
xx xxj xi l zi l 1 1 l 


(23) 


Case III. 


/ [ea(u, x + + fu, *) Su, x + EA(u, x W, W. x + Elw, xx 6w, 

1 j 

.= (f . + F . u i\6u.+ F . - 
\ xi zi — / i y zi 


XX 


dx 


F . — I <5w. + M. 60. 
xi 1 / 1 11 


r 


pu 


+ / p6w + -r— 6u | dx 

r. L R 


u 


(24) 


From the displacement functions, equations (5) and (6) 


= u. y. 


w 


1 ’1,X 

, = w, <p • 

X I r !,x 


( 25 ) 


w, = W. (p. 

XX 1 ^l,xx 
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And the above variations can be written, 


<5u = y. 6u. 

'1 1 

6u, = y. <5u. 

x 'i,x 1 

<5w = (b. 6 w. 

l l 

6w, = 4>. dw. 

X 1,X 1 

dw, = (p. dw. 
xx i,xx 1 


(26) 


Upon substituting equations (25) and (26) into equations (22) through 
(24) and recalling that 0. = w. . the following equations are obtained: 

11 “r Z 


Case I. 


/i EA 

L 


(Vi, + i w k w jVx^,x + kvk,xv) y i,x 6u i 


+ (u.y. w, <p, \<p. dw. 

I j r j,x k k,x / x,x 1 


+ EIw <t> 0. dw.ldx 

j i,xx ],xx if 


- |F . + F . w. . | du. - |F . - F . w. , „ | dw. 

\ xi zi i + 2 J 1 1 zi xi J. + 2 J 1 

+ M i 6w i +2 - 1 [ p h ,5 "’i- pw j‘ (, j,x i 'i 5u iJ d 


Case II. 



( u j r j.x + K w j' , ’k,x*i,x + K U J T k.x r j,xVi . x 5 “ 1 


+/u. y w , <p, ]<p. dw. 

I ] 'j.x k r k,x J i,x 1 


+ EIw. d>. <p. dw. /dx 

J ~1, XX ] ,xx if 


- F 6u - F . dw. - M. <5w. , „ - f p (f>. dw.dx = 0 
xi i zii ii + 2 £ 11 


(27) 


(28) 
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Case III. 


(j EA ( u j Y j,x + i w k w ) ^k, x ^,K + f u kyk,x r j,x^i, x 6u i 
+ (Vj.*"k VxWx 

- (f . + F . -j^Su. -(f 

\ XI zi R I 1 \ 

- f (pep 6w. + U. y. y. 6u. )dx = 0 

L V 1 R Hi i j 


+ EIw . <p . d>. 6w. >dx 

J 1, xx J ) XX 1 / 


w i \ 

. - F . — I <5w. - M. <5w. „ 
zi xx 1 / i ii + 2 


(29) 


For independent virtual displacements, Su. and the equilibrium 
equations are obtained from equations (27) through (29). 

Case I. 

f EA[ u.y. y. + 4-W.W.0, <p . y. + -^-u.u.y 1 y. y. | 

L \ J 3 ’ X 1,x 2 k 3 k » x J» x x » x 2 k J k.x'j.x'i.x J 


+ pw. d>. y. 

J J> x i 


dx-Q.-F.w. _ = 0 
xi zi i + 2 


(30) 


/ 

L 


EAU j r j,x W k' # 'k,x' #, i,x + EI Vi,xx*J.xx - P (*i ) 


dx 


-Q. + F.w. „ = 0 
zi x i i + 2 


(31) 


Case II. 


{ [ EA (Vi.xh.x + K W ) *k, x *i, x \ x + f\ u j'>'k,x->'j,x Y i,x) 


dx 


- Q . = 0 
xi 


( 32 ) 
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■ i ii linn i i 


/ ( EA “j V’lA.x'h.x + EIW j - P *i) dx - Q zi = ° (33) 


Case III. 


/ 


EA (vMV + KVk.Ax Y i,x + fvj\,xV T i.* 


— u. y. y. 

R j 3 1 


dx-Q.-F.-= = 0 
xi zi R 


(34) 


/ ( EAU i ^ ,x w k \,x * i, x + EIW j *1. XX ' *J , xx - P< 4 


Idx 


W) 


- Q . + F . ~ = 0 
zi xi I 


(35) 


Equations (30) through (35) are the final equilibrium equations for the 
beam column element. 

Bifurcation Theory of Instability . Let a solution of the equilibrium 
equations be u®, w®. According to the bifurcation concept of instability, at the 

instability load magnitude there is another set of displacements, arbitrarily 
close to the first set, which also satisfies the equations of equilibrium. Denote 
this second set of nodal displacements by u® + u*, w® + w* . 

Upon substituting this new solution into the equilibrium equations, there 
results for Case I 

/ e 4( u J + u j) n.xn.x + 1(< - w k) ( w j" + w ,‘) *k.**J,X \x 


+ 1 ( u ,® + u, 1 ) ( u® + u? )t, „ % „ % 


2 V k “k/V j J/ k,x j ,x i,x 


+ P 


( w® + w 1 Id). y. >dx 
\ 3 3/ J.x U 


- Q . - F . ( W? 0 + w* ) 
xi ziy _i + 2 _L + ^/ 


= 0 


(36) 
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I 

L 


EA/u? + uA (w° + w 1 )y. 0 0. +El(w° + w 1 '^0 0 

\j ]/ \ k k/ j ,x k,x i,x \] j/ r i,xx j, 


XX 


-P0J 


dx-Q . + F . (w? 

zi xi \ j_+ 2 


+ w . 


• 0 ^ = 0 
I + 2 / 


( 37 ) 


Similar results are obtained for Cases II and HI, 

If these equations are expanded, if the q® state terms (which themselves 

satisfy the equations) are canceled, and if only linear terms in the arbitrarily 
small q* state are retained, the following sets of equations result. 

Case I. 


f EA (ufy. y. + w. 1 w,° 0, 0. y. + 3u,° u?y, y. y. j 

L L V J J,x 1>x J k k,x j ,x i,x k ] 'k.x'j.xh ,xj 


+ pw*0. y. dx - F . Q\ = 0 
j j,x 1 zi 1 


(38) 


I 

L 


EA mx?w° 0. 0 y. +w. 1 u 0 0. 0. y | + EIw*0. 0. 

I ] k i,x k,x j , x j k i,x j , x k,x I j^i.xx^j, 


xx 


dx 


+ f .e\ = o 

xi l 


(39) 


Case II. 


s'- 

L _ 


ea/ u*y. 


V J J 


y: 

X 1,X 


+ wfw° 


k^k.x^j. 


y. + 

X 1,X 




k,x T j ; 


X 


y. \ 

'l,xj 


dx = 0 


(40) 
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/ [EA(^w** liX * kiI Y JiX + w‘u»* t , x * J>x r kiX ) + Etw^i 


dx 


= 0 


(41) 


Case III. 


/ [ EA ( u j\xV + 

w] 

J.x^k.x) +EIW j*t, 


_P 

R 


dx - F . -=• = 0 
zx R 


(42) 


/ 


. <P. 

-xx j ,xx 


dx 


Wi 


+ F . -= = 0 
xi 1 


(43) 


These equations may be expressed in matrix notation for all three cases 
as follows: 


j[K] + [G] + [L p ] + [L F ) {q 1 } = {0} 


(44) 


or 


[K 1 ] {q 1 } ={0} 


where 


[K 1 ] = [K] + [G] + [L p ] + IL f ] 


[K] {q 1 } = 


f EAy. y. 

L 3 ’ X 1>: 


dx 1 


I L 



(45) 
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[G] {q 1 } = 



for all cases. 


The L matrices are different for each case. 


Case I. 


[L p ] M = 


IL F ]{q‘} = 



Case II. 


[L ] = [L J = [0] 
P * 


(47) 


(48) 


(49) 


Case III. 



(50) 
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[L F ]{q‘) = 


Conventional 



Stiffness Matrix K, 

/v 


(51) 


Equation (45) represents the conven- 


tional beam column element stiffness matrix, which is suitable for use when 
there is no interest in nonlinear effects or instability. The elements of this 
matrix may be derived from the general expression, equation (45) . For 
example, 


r.-(‘ -r). *>--('-4 + 4 ) 

from equations (5) and (6). Then, 

1 x x 2 . _ _6_ I2x 

Y l,x = " 17 ’ ^l,x = " 6 L 2 + 6 U * ^l,xx L 2 + L 3 


Therefore, 


K u=/ EAT ltX T 1>x 1 x=f EA (jj) dx =T 


K 33 = 


L L / 6 12x\^ 

/ EI *l,xx'' > l,xx dx = / EI (" 17 + ~L?) dx 


Other elements of K are obtained similarly to yield: 


12EI 

L3 
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[K] = 


u i 

AE 

L 

AE 

L 

0 

0 

0 

0 


u 2 

AE 

L 

AE 

L 

0 

0 

0 

0 


Wi 

w 2 

01 

02 


0 

0 

0 

0 

F xl 

0 

0 

0 

0 

F x2 

12EI 

12EI 

6 El 

6EI 

F 

L 3 

L 3 

“ L 2 

L 2 

zl 

12EI 

12EI 

6EI 

6EI 

F 

-L 3 

L 3 

L 2 

L 2 

z2 

6EI 

6EI 

4EI 

2EI 

Mj 

" L 2 

L 2 

L 

L 

6EI 

6EI 

2 El 

4EI 

m 2 

L 2 

~U 

L 

L 


(52) 


Geometric Stiffness Matrix G. Equation (46) is the general form of the 


geometric stiffness matrix which accounts for the effect of loads existing in the 
element on the stiffness of the element. For example, it is well known that the 
axial load in a beam column has an appreciable effect on the lateral stiffness. 

The geometric stiffness matrix is an adjustment to the conventional stiffness 
matrix to account for such effects. Such matrices have also been referred to in 
the literature as stability coefficient matrices and incremental stiffness matrices. 
Elements of the G matrix will now be derived. 


Evaluation of 


f 


k,x 


4>. y. dx . 
J.x i,x 


For i = 1 , 


y =y -iL^xN _1 
y i,x r i,x dx y L f L 


L 


/ 




x T l,x dx = -L 


I *■ 


k, x 


<p. dx 

J.x 
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when k = 1 , j = 1 



/ 12x 3 

18x 4 36 

x 5 \ 

i 

~12L 3 

18L 4 36 

L 5 

V L4 

L 5 + 5 

L 6 / " 

“ L 

L 4 

L5 + 5 

L 6 J 


0 


12 18 _36_ / 60 90 36 \ _1 6 

~ L 2 L 2 5L 2 “ \ 5 5 5 / L 2 5L 2 ' 

For i = 2 , 




(53) 
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Evaluation of 


<t>. 4>. r, dx 

~i» x j,x k,x 


For i = i 

9 

6x 6x 2 
^i,x *i.x L 2 + L 3 

j = 1, k= i 

f f (p. y A dx = - 77T 
J i,x l,x i,x 51/ 

when j = i, k = 2 , 



(54) 


(55) 
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6 

6 



5L 1 2 

5L 2 



6 

6 


L 

5L 2 

5L 2 


r 1 

f <£, y. dx = 

J 0 r 2,x v k,x r 3,x 

1_ 

1 



10L 

10L 



1 

1 



10L 

10L 



1 

1 



10L 

10L 



1 

1 


L 

10L 

10L 


r 

/ ^3,x\,x r j,x dX 

_ _£ 

_4 

= 


30 

30 


L J 


J_ 




30 

30 



„ 

1 



1 

1 



10L 

10L 



1 

1 


L 

10L 

10L 


i 

f cb <b, y. dx = 

J ^4,x v k,x'j,x 

jl 

_ J_ 

“ 

1 

M 

*©- 

*©- 

0 

30 

30 





_4_ 



30 

30 


L 




Evaluation of j y, y. y, dx 
0 ^ » A J » A 1 * A 

• 



1 1 

Y i,x L ’ ' Y 2,x L 


Then, 


( 56 ) 


( 57 ) 


(58) 


24 



f r, y. y . dx = 

J 0 k »x j,x%x 


{ r k,x 1 'j,x T 2,x dx = 


1 

1 

■ L 2 

L2 

1 

1 

_ L2 

l 2 _ 

1 

1 ~ 

L 2 

~L 2 

1 

1 

_“l 2 

L 2 _ 


[ y k^nj 


[W*] 


(59) 


(60) 


Evaluation of f rf> 0 y dx 
*0 i» x j» x k,x 


f 0 . 0 . y, dx = T0 . 0 

0 1.x j ,x 'k.x 


[*iVk] = 


(61) 


The above matrices are now multiplied by the q® state displacements as 
indicated in equation (46) , k 


W k / < ^k,x 0 j,x Y l,x ClX = [ W i w ° 9 ° e %] 


6 

6 

1 

1 

" 5L 2 

5L 2 

10L 

10L 

6 

6 

1 

1_ 

5L 2 

5L 2 

"lOL 

10L 

1 

1 

__4_ 


10L 

10L 

30 

30 

1 

1 

1 

__4 

10L 

10L 

30 

30 


(62) 
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"k / *k,x*l,x r 2,x & “ [ W ‘ W " °° ] 


6 

6 

1 

1 

5L 2 

5L 2 

" 10L 

10L 

6 

6 

1 

1 

5L 2 

5L 2 

10L 

10L 

1 

1 

_4_ 

1 

10L 

10L 

30 

“30 

1 

1 


_4 

“ 10L 

10L 

“30 

30 


= [w'l 


[*kVj 


(63) 


Similarly, 


w. 


k 


/ ^,x\,x\x dX = [ w!w "‘' ?< ’“] 


6 

6 

5L 2 

5L 2 

6 

6 

5L 2 

“ 5L 2 

1 

1 

10L 

“ iOL 

i 

i 

10L 

“ IOL 


= [w°][0 1 0 k Yj] 


w k / *2,x*k.x 1 'j.x dx " lw ' 1 hM] 


W“ 


/ ■ #> 3,x\,x r i,x dx=[w ° ) [^ 0 k r i] 


(64) 


(65) 


(66) 


W k / V\,x^x d *“ tw ' I (Vk 1 'j] 


(67) 
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< / *i 1 ** J ,x T k.* d3 ‘“ [u * U " 1 


6 

6 

1 

i 

5L 2 

5L 2 

10L 

10L 

6 

6 

1 

1 

5L 2 

5L 2 

"10L 

10L 



= [u°] 

\ *, n 

0 

(68) 


L 




-4 

/ ^2,x^,x\,x dx= [U ° ] [Wk] 


(69) 


L 






Vj 


(70) 


L 




4 

/ *4.x*J,x’'k,x dX=[u * 1 [ 4, 4 

V*] 


(71) 


L 

4 

4 



/ \,x'>'j,x 1 'l,x dx= [U 5 U " 1 

1 

" L 2 

i 

L 2 

= tu°] ^y.rj (72) 


i 

1 




L 2 

L 2 



L 




< 

/ Wux y *.x^ ml ^U* y i y 2l 


(73) 


The above matrices are combined to form the G matrix ( Equation ( 46) ) : 
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3EA[ 


u#, [wi] 


3EA[u°] |^7 k y^ y 2 J 
EA[w°] 


[G] = 


EA[w°] ^ 2 \ r j] 


EA[w°] ^ 3 <^ k y. 


EA[w°] [^Vj] 


w| w| 0i 02 

EA[w°]^ k ^ y i J 


EA[w°][^ k 0.y 2 J 

EA[u°] ^i^jT k J 


[u#1 [ VjV 


j EA[u»][^0.T k J 

| EAlu 0 ]^ * 4 *jTj. j 


Load Behavior Matrix. The effect of applied load behavior on the element 
stiffness is obtained by adjusting the K matrix with the L matrix, where 

[L] = [L p ] + [LjjJ . (75) 


The L matrix will now be derived from equations (47) through (51) for each of 
the three load behavior cases. 


Case I. 


5 pw j^j x r i dx= / p^-l) w K~iJ + l0 +0 K" 1 + 4 l“ 3 l*) 

A J J > T . 


i/ 6x 6x“ 

Hl 2 " l3 


) + d(fr ■ d * 


= - |pw{ + |pw 2 * - ^P L£, 1 + ^P l0 2 • 
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Then, 


1 1 1 i 

/ PW !<P. 7 2 dx = -pw 1 + -pw 2 +— pL0j - ^ pL0 2 


0 

0 


0 

0 


1 i _L_ 

2 2 " 12 12 

1 i L_ _L_ 

2 2 12 12 


L 

P 


P 


0 

0 


0 

0 


0 

0 


0 

0 


0 

0 


0 

0 


0 0 


0 0 0 0 


0 


0 0 0 0 0 


(76> 


L may be written directly from equation (48) in view of equations (38) 


and (39) ; 


[l f ] = p 


0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 


0 

0 

0 

0. 

0 

0 


-F 

0 

F 

0 

0 

0 


zl 


xl 


0 

-F 

0 

F 

0 

0 


z2 


x2 


(77) 


where F's are proportionality constants between the indicated forces and p. 

The magnitudes of the applied forces are assumed to remain in a constant ratio 
to each other and to the lateral load, p, during loading of the element. This 
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assumption is made here and in the Stability Criterion section for illustrative 
purposes. The assumption is easily altered to study buckling under pressure 
with specified nodal forces or buckling under nodal forces with specified 
pressure. 

Case II. 


[L ] = [L ] = [0] 
p F 


(78) 


Case III. 


/ I [ u K‘-i) +u Ki) 




dx 


JL r 

R • l 




_ X 
R 


( x -r + ii?- 1 ui ns -5i? 


x- ) „} + f _ _5l) u l 

j( L -L + i) u l + (i-i) 4_ 


= X 

R 


f u*y 
J T R j 7 3 


, pL 

: dx = T 


1 i 1 j 

Le + 3^1J 


[L ] = p 
P 


_L 

3R 

x 

6R 

0 

0 

0 

0 


_L_ 

6R 

_L 

3R 

0 

0 

0 

0 


0 

0 

0 

0 


0 

0 


0 0 

0 0 

0 0 
0 0 


0 0 0 0 
0 0 0 0 


(79) 
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Lj, may be written directly from equation (51) ; 



31 



Coordinate Transformation 


The developments thus far have been in a coordinate system which was 
oriented so that the x-axis coincides with the element centerline. To combine 
several elements for solution of a particular problem, it is necessary to obtain 
the stiffness matrices of the elements in a common or global coordinate system. 
This system will be denoted by x, z as shown in Figure 5. 



FIGURE 5. COORDINATE TRANSFORMATION 

The figure shows that vectors in the two coordinate systems transform 
according to the relationships 

m = m cos P + h sin P 

n = n cos P - m sin p 

CO = CO . 
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Writing in matrix notation, we see that 



{m} = [T] {m} . ( 81) 

Thus the displacements and forces of the element transform according to 
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( 82 ) 


(83) 


If the nodal forces and displacements are related by the stiffness K 1 
( equation (44 ) ) , then 



(84) 


where K 1 has been rearranged to conform to the Q and q matrices. Upon sub- 
stituting equations (82) and (83) into equation (84) , 


T 1 0 

( Qi f 

T 1 0 

f 

1 

1- _ 

| 

= [K 1 ] 

_0 ! T_ 

f Q2 ) 

_° 1 T 




T i 0 

-i 

T 1 0 

_0 | T_ 

[K 1 ] 

L _ 

_° ! t _ 



But in this case, 




-1 

I — 

T j 0 


T , 0 

_0 ! T_ 


0 | T_ 


( 85 ) 
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I 

i 0 



[K 1 ] 




( 86 ) 


( 87 ) 


where K is the stiffness matrix of the element in the global coordinate system. 


Master Stiffness Matrix 


The stiffness matrix for complete structure is obtained by combining the 
element stiffness matrices in the global coordinate system. This section 
describes the method of accomplishing this combination. 

Consider part of an overall structure schematically represented by 
Figure 6. The Q.'s in Figure 6 are intended to represent the total load vector 



FIGURE 6. GENERAL STRUCTURE 
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at the point of application, 


Q i= F ,i 


Load-displacement relationships for each of the three elements in the 
global coordinate system are expressed from equation ( 87) as follows: 


A 1 -- 2 ! 

Q . J 

ai- i 


I _ 

TC I K 

ai - 2, i-2 , ai-2, i-1 

F 

K ai- 1, i-2 J K ai- 1, i- 1 


Vi 


bi-i/ = ^i-l, i-i | K bi-i, i 

bi \ *Sbi, i-i ! ^bi» i 


Vl 


Q . , 1 

<3i + i 


1 K . ... 

I Cl, 1+ i 


^ei+i, i I ^ci+i, i+1 

i 


■ + i\ . 


By combining the above relationships with the fact that q. is the same for 
both elements b and c, and noting that 


Q i = Q bi + Q cl 


and that similar relationships hold for the other nodes, the master stiffness 
matrix is obtained: 
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0 


0 


0 



— 

Vi.i-i+Vl. 

1 

£ 

1 

0 

Vi-1 

V..+K . 

bn cii 

0 

0 

K , 4 , 

Cl + 1 , l 

0 

0 

0 



(89) 

Note that this is a relationship between the externally applied loads 
and the nodal displacements of the -assembled structure in the global coordinate 
system. 

Stability Criterion 

Elements with stiffness matrices of the type given by equation (44 ) may 
be assembled into a master stiffness matrix to represent a structure subjected 
to the critical (buckling; magnitude of applied loads. Thus, 


[K] + p [K 0 ] 


{qi}={0} 


where. 


(90) 


p[K 0 ] = [G] + [L ] + [L p ] , (91) 

and the null matrix of applied external loads indicates that the structural stiffness 
has vanished under the critical load magnitude (a physical interpretation of 
instability) . Boundary conditions may be applied to reduce the size of the set of 
equations (90) and the new reduced set again denoted by the equations (90). 

The K, G, and 17 matrices above are obtained by assembly of element 

~ ~ ~p 

matrices as described in the preceding sections. The L^, matrix is more 

conveniently obtained by direct application of equations (48) and (51) to each 
node of the assembled structure in the global coordinate system. 

A nontrivial solution of equations ( 90) will exist only when the determinant 
of the matrix K + p K 0 vanishes ; 


K + pK„ 


= 0 . 


(92) 


This is an eigenvalue problem where the magnitudes of applied load, p, at which 
instability will occur are the eigenvalues. 
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Primary Equilibrium State 


It should be noted at this point that the foregoing solution is contingent 
upon a knowledge of the primary equilibrium state. That is , the G matrix can- 
not be formulated until the q° displacements are known in terms of the applied 
loads. Given a set of applied loads, the nonlinear equilibrium equations (30) 
through (35) may be solved for the q° displacements. An iterative procedure 
for incorporating this solution into the stability problem is given in Reference 14. 

Another approach which is consistent with many classical stability prob- 
lem solutions is to linearize equations (30) through (35) for the purpose of 
obtaining the primary state solution. If this is done these equations become for 
all three cases 


/ 


EAu?y y. dx - Q . = 0 

J J , X 1, X XI 


(93) 


/ 




EIw ?0. 

J i,xx 


0. 

J.xx 


- P0. 


dx - Q . = 0 
zi 


(94) 


All terms in the above equations have previously been evaluated and 
expressed in matrix form except the p term which is as follows: 



1 

f p<p 2 dx = ~pL- 
0 2 


L 

f pcp 3 dx = 

0 



pl-x + 





dx 
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JL- 1 

/ p<£ 4 dx= — pL 2 L 


or 


{p<£} = p 


0 

0 

L 

2 

L 

2 

it 

12 

12 


(95) 


Equations (93 , 94) may then be expressed in the form 


[K] {q 0 } ={p 4 >) +{Q} . 


(96) 


From this expression for each element, the master stiffness matrix for 
the entire structure is assembled as before. Then for the complete structure, 


{q 0 } = [K] 1 {p </> + Q} 


(97) 


A further simplification for obtaining c^° which has also been used 
extensively in classical problems is the use of a membrane solution for q°. The 
structure is assumed to take no loads in bending before buckling. For many 
practical problems such a solution can readily be obtained by inspection, ele- 
mentary equilibrium considerations, or from the literature. For simplicity this 
approach is used in the example problem of this report. 


APPLICATION OF THE THEORY 


This section is intended to give a step by step approach for applying the 
theory to the solution of practical problems. For problems with several elements 
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either a general computer program would need to be developed or a relatively 
simple computer program for each application could be developed as was done 
for the example to follow later. 

One begins by dividing the structure to be analyzed into discrete elements, 
the number depending upon the desired accuracy. From the known properties 
and loading on each element the K matrix can be calculated from equation (52) 
and the L p and L^, matrices determined from the appropriate equations (76), (77), 

(78), (79), or (80). Note that the critical magnitude, p, of the applied loading 
remains as the unknown to be determined. 

The q° primary state is now determined from a known membrane solution 

or from equation (97). The K matrix in equation (97) is assembled by using 

equations (81) and (86) to transform elements to the global coordinate system 

and equation (89) for combining the elements. Boundary conditions are applied 

and equation ( 97 ) is solved for q°. Once the q° state is determined the G 

/■>*> 

matrix can be found from equation (74), using the definitions (53) through (61). 

The K 1 matrix defined by equation (44) is now known for each element. 
These are assembled into a master stiffness matrix by equations (86) and (89) , 
and boundary conditions are applied to the resulting set of equations to obtain 
equations (90). Eigenvalues of the characteristic determinant, equations (92), 
which is obtained from equation (90), are the magnitudes of the buckling load. 

The mode shape for each buckling load can be determined by substituting each 
eigenvalue in turn into equations (90) and solving for the relative amplitudes of 
displacements, qL. 


EXAMPLE PROBLEM 
Circular Arch With Uniform Pressure 


The above application procedure was applied to the uniform circular arch 
shown in Figure 7. The known membrane solution for a circular arch under 
uniform pressure is 
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( 98 ) 


[q°] 


= [u°w°] - [u$ 

_ jdR [~,L _ Ij 

AE L 2 "2 


u| Wj 
- R 


w 2 ° e° t e§] 

-Roo] 


for all elements. 

The G matrix is formulated as indicated by equation (74) : 


3EA [u°] 


Wi 




1 

1 

L 2 

L 2 

1 

1 

L 2 

" l 2 _ 


[- 1 o 

3EA[u»][ Tkrj y 2 ]=3EK ["! 


EA [u°] 





6 

6 

1 

1 

5L 2 

5L 2 

10L 

10L 

6 

6 

1 

1 

_ 5L 2 

5L 2 

~ 10L 

10L 


= pR 


_6_ _6_ _1_ _i_ 

5L 5L 10 10 
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RADIUS = R = 100 inches (2.54m) ^ 

AREA = A = 0.628318 in 2 (4.05 cm 2 ) 

MOMENT OF INERTIA = I = 0.314159 in* (13.1 cm 4 ) 

YOUNG’S MODULUS = E = 10 7 ps i (6.9 x 10 10 N/m 2 ) 

FIGURE 7. UNIFORM CIRCULAR ARCH 


Other elements of the G matrix are calculated similarly to yield for all 
elements 


[G] — pR 


u i 

u 2 

Wi 

w 2 

01 

02 


_3 
" L 

3_ 

L 

0 

•o 

0 

0 

F xl 

3_ 

L 

3_ 
" L 

0 

0 

0 

0 

F x2 



_6__ 

_6_ 


1 


0 

0 

" 5L 

5L 

10 

10 

zl 



6 

_6_ 

J_ 

1 


0 

0 

5L 

“ 5L 

■ 10 

“10 

F z2 



i 

JL 

4L 

_L_ 

Mj 

0 

0 

10 

" 10 

" 30 

30 



1 

J_ 

_L 

4L 

m 2 

0 

0 

10 

10 

30 

~ 30 


(99) 


where the four terms in the upper left corner can be traced directly to the 
u J term in equation (11). 



Now for all elements the K matrix is given by equation (52) , the L 

P 

matrices for the three load cases are given by equations (76), (78), and (79), 
respectively, and all L matrices are null. 

rs ' / r 

The K, G, and L matrices are now all rearranged to the order, 

r^ r^f Jr 


[K 1 ] = 


u i w i d l u 2 w 2 9 Z 

I ”” 


K fi 


K|i 




K 22 


xl 


zi 


I 




J 


x2 


z2 


M_ 
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All element matrices are now rotated to the global coordinate system, 
x, y by the operation, 


m 


l— — ■ 

I 

T 0 

[K 1 ] 

T 0 

T 


1 — 

o 

H 

1 

1 

o 

H 

1 


IK] 


( 101 ) 


as indicated by equation (86) . The T matrix is of course calculated separately 

for each element as indicated by equation (81) where (3 is the angle measured 
clockwise from the element x axis to the x axis. 

The element matrices are now combined to form the master stiffness 
matrix as indicated by equation (89) . The size of this matrix is reduced by 
applying the boundary conditions , 

u 1= Wi = e ± = ^ = Wj = 0 = 0 , (102) 

where u^, w^, 6 ^ are displacements of the right end of the last element. From 
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this reduced matrix the characteristic determinant (92) is formed and eigenvalues 
are calculated. The above operations were performed on a computer for idealiza- 
tions of 2, 3, 6, 9, and 12 elements. A hand calculation was also made for the 
two-element arch. Results of the analysis are shown in Table I. As can be seen 
from the table, comparison with known results is excellent and there is rapid 
convergence to the exact solutions. The exact solutions were obtained from the 
analysis by Wempner and Kesti; [24]. The 12-element solution for Case II has 
also been previously obtained by the finite element method in Reference 14. 


TABLE I. COMPARISON OF ARCH BUCKLING LOADS 


Number 

of 

Elements 

Buckling Pressure (lb/ in) 

Case I 

Case II 

Case III 

2 

94. 1 

94. 1 

94. 1 

3 

57. 1 

64. 7 

67.4 

6 

57.4 

62. 6 

64.4 

9 

56.6 

61. 9 

63. 6 

12 

56.4 

61. 6 

63.3 

Exact solution 

56. 87 

60. 95 

63.46 


CONCLUSIONS 


The results of the example problem above show that the theory 
developed here may be quite useful in solving practical engineering stability 
problems. 
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Ill 


The implicit assumption of uniform shear on the element mentioned 

in Development of Element Stiffness Matrix section does not appear to adversely 

affect the results. However, if desired, this assumption could possibly be 

eliminated by the technique introduced by Pian [19]. Also, the u 2 term of 

, x 

equation ( 11) apparently has a negligible effect on the results for the particular 
example solved, since the exact solutions used for comparison do not contain 
effects of this term. 

A very important item to note is that the present theory applies to 
unconservative systems (represented by Case I) as well as conservative 
systems since the principle of virtual work was used. It is known, however, 
that certain unconservative systems are stable in the static sense (considered 
here), but unstable in the dynamic sense. (See Reference 25, page 152 
through 156. ) Thus the present analysis could not be expected to adequately 
treat this type of problem, and great care should be used in applications to 
unconservative systems. 

Fruitful extensions of the present theory would likely be found in the 
area of curved and three-dimensional elements and in a general automation of 
the theory application. 


George C. Marshall Space Flight Center 

_ National Aeronautics and Space Administration 

Marshall Space Flight Center, Alabama 35812, May 1, 1968 
981-10-10-0000-50-00-008 
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